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A review of the most popular Linear Multistep (LM) Methods for solving 
Ordinary Differential Equations numerically is presented. These methods are 
first derived from first principles, and are discussed in terms of their order, con- 
sistency, and various types of stability. Particular varieties of stability that may 
not be familiar, are briefly defined first. The methods that are included are the 
Adams-Bashforth Methods, Adams-Moulton Methods, and Backwards Differ- 
entiation Formulas. Advantages and disadvantages of these methods are also 
described. Not much prior knowledge of numerical methods or ordinary differ- 
ential equations is required, although knowledge of basic topics from calculus 
is assumed. 

1 Linear Multistep Methods 

As opposed to one-step methods, which only utilize one previous value of the 
numerical solution to approximate the subsequent value, multistep methods approx- 
imate numerical values of the solution by referring to more than one previous value. 
Accordingly, multistep methods may often achieve greater accuracy than one-step 
methods that use the same number of function evaluations, since they utilize more 
information about the known portion of the solution than one-step methods do. 

A special category of multistep methods are the linear multi-step methods, where 
the numerical solution to the ODE at a specific location is expressed as a linear 
combination of the numerical solution's values and the function's values at previous 
points. For the standard system of ODEs, y' = f(t,y), a linear multistep method 
with k-steps would have the form: 
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where otj, f3j are constants, y n is the numerical solution at t — t n , and f n = f(t n , y n ). 
For the rest of this discussion, we will make the assumption that / is differentiable as 
many times as needed, and we will consider the scalar ODE y' = f(t, y) for simplicity 
in notation. The generalization to systems of ODEs is straightforward. 

It is important to note that in the above expression, all of the previous integration 
steps are assumed to be equally spaced, although it is possible to generalize these 
schemes to have variable step-sizes. Also, note that if /?o = , the scheme is explicit 
(because it does not depend on f n ), and otherwise the scheme is implicit. We are now 
ready to examine some of the most popular families of linear multistep methods. 

1.1 The Adams Family 

1.1.1 The Adams-Bashforth [AB] Methods (Explicit Adams Methods) 

The most widely utilized linear multistep methods used for nonstiff problems 
are the Adams-Bashforth methods, which are members of the Adams Family that 
are explicit. The spirit of the Adams-Bashforth technique is rooted in the Stone- 
Weierstrass Theorem. 

Theorem 1 - Stone- Weierstrass Theorem - Let f(t) : H. — > C be continuous on 
t E [a,b]. For all e > 0, 3 a polynomial 4>(t) 3 \ \f(t) - <f>(t)\\ < e. 

In other words, any continuous function can be approximated to an arbitrary accuracy 
by a polynomial; generally, the more demanding the accuracy of the approximation, 
the higher the order needed of such a polynomial. 

With the Stone- Weierstrass Theorem in mind, we start with the ODE in question: 
y' = f(t, y), and we integrate both sides to obtain: 



If we could integrate f(t,y(t)) analytically, we (likely) would not need to resort 
to numerical methods to determine the solution to the ODE. If we cannot integrate 
f(t,y(t)) analytically, according to the Stone- Weierstrass Theorem above, we can 
approximate it with arbitrary accuracy by a polynomial 4>(t), and since all polynomi- 
als can be integrated analytically, we have an obtainable, fair approximation of the 
solution to the ODE: 
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Now to ensure that our approximation is reasonable, we insist that 0(t r 
/(t n _i) for a reasonable number of integer values i. For example, setting 0(t„_i) to 
be the constant /(t n -i) will result in the scheme: 

y{t n ) « y(i n _ 1 ) + /i/ n _! (4) 

Note that this scheme, which is also known as the 1-step Adams-Bashforth Method, 
is simply the classic Forward Euler (FE) method. In terms of equation ([T]), this scheme 
has «i = —1, /3 = 0, Pi — 1 and Pj, ctj = for j > 1. 

Let us now construct the 2-step Adams-Bashforth scheme. We first need an inter- 
polation polynomial <p(t) such that 0(t„,_j) = f{t n _i) for i — 1,2. The desired linear 
function is displayed below: 

f(t, y) « 0(f) = /(f n _ 2 ) + /(t ?- l} ~ f (tw ' 2) (t - f n _ 2 ) (5) 
Together with ([3]), we have that 

f{t n -l) ~ f{tn-2) {t — tn-2) 21 



V(tn) ~ V{tn-l) + 



f(t n - 2 )t + 



t n -l ~ t n -2 



(6) 



tn-l 



y(t n ) w y(i n _0 + /i ( - J/(tn- 2 ) ) (7) 



So in accordance with the line above, we can define our 2-step Adams-Bashforth 
scheme to be: 

3 1 

Vn = Vn-1 + H^fn-l - 2^- 2 )- ( 8 ) 

This can also be expressed in terms of (JTJ with a,\ = — 1, /3q = 0, (3\ = §, /3 2 = — \ and 
ctj = 0, for all i > 1. 

Continuing in this manner, we can construct fc-step Adams-Bashforth methods 
by interpolating / through k previous points: t = t„_i, t n _ 2 , . . . £ n -fc- Such a scheme 
could be derived by constructing a degree < k — 1 polynomial such that 0(t n _ i ) = 
f(t n -i) for z = 1,2,... A;, and integrating it as in (J3J) , then replacing y(t n ),y(t n ^i) 
and f{t n -i) with y n , y n _i and / n _j respectively. As shown below, the resultant /c-step 
Adams-Bashforth method can be expressed in the form of (0Q) with a\ = —1, (3q = 0, (3j 
defined as displayed below for 1 < j < k, and cxj = 0, Pj+k-i = for j > 1: 

k 

Vn = Vn-1 + h ^2 Pifn-j, (9) 
3=1 

where 

fe-i / . \ ,1 / 
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It is important to mention that for such schemes, k starting values must be given. 
If only the initial condition is provided, the other k — 1 points can be determined by 
a different scheme (for example, a Runge-Kutta method of the same order). Adams- 
Bashforth methods also tend to have small regions of absolute stability (to be dis- 
cussed later), and this inspired the construction of implicit Adams methods (called 
Adams-Moulton methods) which are the topic of the following discussion. 

1.1.2 The Adams-Moulton [AM] Methods (Implicit Adams Methods) 

The difference between Adams-Moulton and Adams-Bashforth methods is that 
Adams-Moulton methods use an interpolating polynomial of degree < k rather than 
< k — 1 , and it includes / at the unknown value t n as well. A fc-step Adams-Moulton 
scheme can be expressed in the form of ([!]) as follows: 



It is apparent that when k = 1 and /3\ = we have the classic Backward Euler 
(BE) method. Likewise, if k = 1 and fix ^ we have the implicit trapezoidal method. 

Adams-Moulton methods have smaller error constants, use less steps, and have 
larger stability regions than their Adams-Bashforth counterparts (of the same or- 
der). However, AM methods using more than one step tend to have smaller regions 
of absolute stability than other implicit methods such as Runge-Kutta methods (in 
fact, they tend to be bounded, which often defeats the purpose of using an implicit 
scheme). Adams-Moulton methods have smaller error constants, use less steps, and 
have larger stability regions than their Adams-Bashforth counterparts (of the same 
order). However, AM methods using more than one step tend to have smaller regions 
of absolute stability than other implicit methods such as Runge-Kutta methods (they 
tend to be bounded, which often defeats the purpose of using an implicit scheme). 
Adams-Moulton methods have smaller error constants, use less steps, and have larger 
stability regions than their Adams-Bashforth counterparts (of the same order). How- 
ever, AM methods using more than one step tend to have smaller regions of absolute 
stability than other implicit methods such as Runge-Kutta methods (they tend to be 
bounded, which often defeats the purpose of using an implicit scheme). 

An alternative family of implicit linear multistep methods is the family of Back- 
wards Differentiation Formulas, which are the topic of the next section. Such schemes 
are in fact the most popular methods for stiff problems. 

2 Backwards Differentiation Formulas [BDFs] 

In contrast to the linear multistep schemes in the Adams Family, who are derived 
by integrating an interpolating polynomial (p(t) that approximates /, the BDF meth- 
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ods are derived by differentiating an interpolating polynomial <p(t) that approximates 
y (one such that ip{t n _i) = y{t n _i) for i — 0, 1, 2, . . . k), and setting the derivative at 
t n to be equal to f(t n ,y n ). 

For example, the 1-step BDF method is derived as follows. We first construct 
the interpolating polynomial tp(t) that approximates y, with ip(t n ^i) = t/(£ n -i) for 
i = 0,l. 

y(t) « <p{t) = y{t n ) + (t- tn ) V ^-_f n - l) (12) 

tn tn—1 

Upon differentiation, we get: 

y'(t) = f(t, y) » y/(f) = ^ ~ ^) , (13) 

£71 "n— 1 



We can then use the approximation in ( 1131) as inspiration to construct our 1-step BDF 
method, by setting (p'(tn) = f(t n ,y n )- 

y(K) -»(i_J = 

This is in fact, not surprisingly, equivalent to the Backward Euler method when 
rearranged. 

Similarly, we can construct a fc-step BDF by generating the fc-degree interpolating 
polynomial: 



y(t) « tp(t) = y n + \{t - t n )Vy n + -^{t - t„)(t - t n ^)V 2 y n + ... + p^(t - t„) . . . (t - t n _ k+1 )V k , (15) 

where V* is the backward difference operator: 

VV = y n (16) 

V l y n = V l ~V- V^Vn-i- (17) 
Then upon differentiating, and setting (p'(t n ) = f(t n ,y n ) we get: 

Y l -V l y n = hf(t n ,y n ), (18) 

which can be transformed to match the general expression of (CQ), with f3j = for 
j > (note that this makes them implicit schemes): 

Vn = -}, ka { y n -i + /i/3b/„. (19) 

i=l 
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It is noteworthy that the backward differences J \7 l y of y approximate the true 
derivatives of y (i.e. V*y ~ h k y^>). 

What makes BDF methods powerful is their unique and convenient region of 
absolute stability giving rise to its L-stability (to be discussed later), which is of 
particular importance for stiff problems. It is for this same reason that BDF methods 
with k > 6 are not used, as these methods have a region of absolute stability that 
crosses the negative real axis, thereby disqualifying them from being classified as 
A-stable. 

3 Order and consistency of Linear Multistep Meth- 
ods 

Investigating convergence of linear multistep methods is quite different from that 
of non-linear one-step methods (such as the Runge-Kutta) methods. For Runge-Kutta 
methods, O-stability is automatic, and investigating the order can be cumbersome. 
For linear multistep methods, O-stability is not necessarily automatic, and needs to be 
confirmed for each scheme. Contrarily, and again unlike the Runge-Kutta methods, 
investigating the order of linear multistep methods is rather straightforward. We will 
start this section by investigating the order of linear multistep methods. 

3.1 Order 

To begin, let us define the linear operator Jzf^ [?/(£)], where y{t) is an arbitrarily con- 
tinuously differentiable funciton on [0,b]: 

k 

^h[y(t)) = Y^[a jy (t - jh) - hfyy'it - jh)) (20) 

3=0 

This expression is based on Eq. (pQ). Recalling that y' = f{t,y(t)), we can write the 
above expression in the following way: 

k 

■&h[y(t)] = $> j2 /(t - jh) - h/3jf(t - jh, y(t - jh))} (21) 

j=0 

which becomes, after expanding y(t — jh) and f(t — jh, y(t — jh)) in a Taylor series 
about t and simplifying: 

S? h [y{t)] = C y(t) + dhy'it) + ... + C q h 9 h"(t) + ..., (22) 

where, 
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Co = aj , and 

3=0 



C; 



1,2,3,... 



(23) 
(24) 



. 3=1 V ' 3=0 

Now, the order of the method is p if the local truncation (or discretization) error 
is d n = 0(h p ), which is given by: 



d„ 



h 



(25) 



where y{t n ) is the exact solution at t — t n . So using the version of the expression 
given by Eq. (1221) . we get that the method is order p if: 

C = C 1 = ... = C p = 0,C p+1 ^0. (26) 
Combining this result with (122]) . we get that 

d n = C p+1 h p y^ +1 \t n ) + 0(h? +1 ), (27) 

where C p+ i is the error constant of the scheme. 

It can then be shown that Adams-Bashforth and BDF methods are of order k 
(where k is the number of steps) , while Adams-Moulton methods are of order k + 1 
(with the exception of the case where the scheme is completed in a single step with 
f3\ = 0, as in Backward Euler, which is order k — 1). 

Calculating each value of C q in ( 1221) can be tedious though. The easiest way to 
check for consistency (p > 1) is to use the fact that a method is consistent if and only 
if 

k k k 

<*3 = and ^3 + J2 & = °" ( 28 ) 

3=0 j=l 3=0 

This result can also be demonstrated in terms of the characteristic polynomials of 
the recurrence relations arising from the expression for the numerical scheme: 



3=0 



3=0 

In particular, the scheme is consistent if and only if p(l) = 0, p'{l) = cr(l)- 



(29) 
(30) 
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4 Stability of Linear Multistep Methods 

4.1 O-Stability 

Here we are looking at Eq. ([1]) in the limit as h — > 0, rendering this equation 
into the form: 

a k y n _ k + a k _iy n ^ k+ i + ... a y n = 0. (31) 

This recurrence relation has the characteristic polynomial defined earlier. Due 
to consistency (see previous section), we know that £ = 1 is a root of Eq. (I3TI) . If the 
rest of the roots of the equation are distinct, then we have a solution of the form: 

fc-i 

!fe = 5> + CP + <*(!)"■ (32) 
i=i 

If £i = £ 2 is a double root, then the solution will have the form: 

fc-i 

Vn = « + # + C ^ + + C 2<2 • (33) 
i=3 

Similarly, if £i = £2 = £3 is a triple root, then the solution will have the form: 
fc-l 

Vn = Yl C * + & + C ^ l T + + C 2<2 + C3Tl(n - (34) 

i=4 

We can see that if |£«| > 1 , then our solution will diverge as n gets large. Likewise, 
if = 1 is not a simple root, we will again have divergence. Therefore, a linear 
multistep method is 0-stable if and only if all roots of the equation = satisfy 
l&l < 1 , where if = 1, then ^ is a simple root, for 1 < % < k. Now by a theorem 
sometimes referred to as the Dahlquist Theorem, if this root condition is satisfied, 
and the method is accurate to order p, and the initial values are accurate to order p, 
then the method is convergent to the order p. 

We can further specify the strength of the stability of the scheme by defining 
strongly stable as meaning all roots of = have the property |£| < 1 with the 
exception of the one root £j which equals 1. A scheme can then be defined as weakly 
stable if it is 0-stable, but not strongly stable. 

4.2 Absolute Stability, A-stability and L-stability 

By applying Eq. (p^) to the test equation y' = Xy and letting y n = £ n , we find 
that £ must satisfy p(£) — h\a(£). Now to address absolute stability, we define the 
stability polynomial as ( = p(£) — h\cr(£) and note that the absolute stability region 
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is the region of values of Xh such that \y n \ does not grow with increasing values of n. 
For this condition to be met, all roots £j of £ = must satisfy \^\ < 1. 

As for A-stability, recall that a numerical scheme is A-stable if its absolute stability 
region contains the entire left half of the complex plane (i.e. it contains Ke(Xh < 
0)). It turns out that explicit linear multistep methods cannot be A-stable, and 
that A-stable linear multistep methods with order greater than 2 do not exist. The 
most accurate LM method (method with smallest error constant) is the second-order 
implicit trapezoidal method (with error constant C3 = i). We conclude that A- 
stability is rare in linear multi-step methods. 

After applying our numerical scheme to the test equation y' = Xy, we can rearrange 
our expression to obtain an expression of the form y n = R(z)y n -i, where z = Xh. If 
R(z) — > as Re(z) — > —00, we say that the scheme is L-stable (or has stiff decay. 
This explains why the Backward Euler method works better than the Trapezoidal 
Rule for some problems, even though the Trapezoidal Rule is of higher order! 



10 



